Path independent integrals to identify localized plastic events in two dimensions 
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\ We use a power expansion representation of plane elasticity complex potentials due to Kolossov 

and Muskhelishvili, to compute the elastic fields induced by a localized plastic deformation event. 
' Far from its center, the dominant contributions correspond to first order singularities of quadrupolar 

and dipolar symmetry which can be associated respectively to pure deviatoric and pure volumetric 
plastic strain of an equivalent circular inclusion. Constructing holomorphic functions from the 
' displacement field and its derivatives, it is possible to define path independent Cauchy integrals 

which capture the amplitudes of these singularities. Analytical expressions and numerical tests 
on simple finite element data are presented. The development of such numerical tools is of direct 
interest for the identification of local structural reorganizations which are believed to be the key 
(JJ) ' mechanisms for plasticity of amorphous materials. 
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I. INTRODUCTION 

.. 

Plasticity of amorphous materials has motivated an increasing amount of studies in recent years. In the absence 
I ' of underlying crystalline lattice in materials such as foams, suspensions or structural glasses, it isgenerally accepted 
y [ that plastic deformation results from a succession of localized structural reorganizations [ll H, [J Such changes of 
local structure release part of the elastic strain to reach a more favorable conformation and induce long range elastic 
Q I fields. The details of such local rearrangements and of the internal stress they induce obviously depend on the precise 
structure of the material under study, and its local configuration . However, the important observation is that outside 
the zone of reorganization, a linear elastic behavior prevails. Therefore, elastic stresses can be decomposed onto a 
multipolar basis and independently of the material details, it is possible to extract singular, scale-free, dominant terms 
■ which can be associated to a global pure deviatoric or pure volumetric local transformations of an equivalent circular 
J — ' inclusion. In particular, the elastic shear stress induced by a localized plastic shear exhibits a quadrupolar symmetry. 
Cn| ' This observation has motivated the development of statistical models of amorphous plasticity at mesoscopic scale 
' based upon the interaction of disorder and long range elastic interactions!^ S S ^^^^ same spirit statistical 

. models were also recently developed to describe the plasticity of poly-crystalline materials [H. Several numerical 
' studies have been performed recently to identify these elementary localized plastic events in athermal or molecular 
OO , dynamics simulations of model amorphous material under shear [lol . ITl| . 

' The question remains how to identify and analyze these transformation zones. In analogy with the path independent 
Rice J- integral [l3 developed to estimate the stress intensity factor associated to a crack tip stress singularity, we aim 
here at capturing the stress singularity induced by the local plastic transformation which can be treated as an 
Eshelby inclusion [l^. In two dimensions, we develop a simple approach based upon the Kolossov-Muskhelishvili 
^ . (K&M) formalism of plane elasticity This is an appealing pathway to the solution since these zones will appear 
as poles for the potentials, and hence Cauchy integrals may easily lead to contour integral formulation which are 
independent of the precise contour geometry, but rather rely on its topology with respect to the different poles which 
are present. 

Although these techniques have been mostly used in the context of numerical simulations in order to estimate stress 
intensity factors from finite element simulations, they are now called for to estimate stress intensity factors from 
experimentally measured displacement field from e.g . digital image correlation techniques. In this case, interaction 
integral techniques or least square regression (la| techniques have been applied. Noise robust variants have also 
been proposed [17|- These routes could also be followed in the present case. 

Though the present work is restricted to two dimensions due to the complex potential formulation, similar questions 
can be addressed for the three dimensional version of this problem using the same strategy but a different methodology. 
In the following, we briefly recall the K&M formalism, and we give analytic expressions of contour integrals allowing 
to capture the singular elastic fields and we present a few numerical results based on a finite element simulation 
supporting our analytical developments. 
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In Section 2, we present the theoretical basis of our approach in terms of singular clastic fields, while Section 
3 introduce the contour integral formulation. In section 4, a numerical implementation based on finite element 
simulations is presented together with the results of the present approach. This application allows to evaluate the 
performance and limitations of the contour integral procedure and check the detrimental effect of discreteness. Section 
5 presents the main conclusions of our study. 



II. POTENTIAL FORMULATION 



In two dimensions the Kolossov and Muskhclishvili potentials (K&M) can be used to write the elastic stress and 
displacement fields U and cril4|]. Using a complex formulation, we introduce the elastic displacement \J = Ux + iUy 
and the stress tensor field through two functions, the real trace 5*0 = <Jxx + fyy and the complex function S = 
ayy — Gxx + lioxy- In the framework of linear elasticity, balance and compatibility equations can be rewritten as: 

5o,. - S,- = , (1) 
^0,.- = , (2) 

where z — x ^ iy is the complex coordinate and the notation A,x is used to represent the partial derivative of field 
A with respect to coordinate x. Note that we assumed zero surface density force and that equation ^ is here the 
classical Beltrami equation which expresses the kinematic compatibility condition in terms of stress. The general 
solution to these equations can be obtained through the introduction of two holomorphic functions Lp and ill, called 
the K&M potentials. The displacement and the stress field can be written [iJl 
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Kip(z) — Z(p'{z) — ^j{z) 



where fi is the elastic shear modulus and k 
being the Poisson's ratio. 



(^'(z) + (^'(z) 



2[V(z) + ^'(z)] 



(3) 
(4) 
(5) 



(3 — 4:v) for plane strain and k = (3 — + for plane stress, v 



III. PLASTIC INCLUSION AND SINGULARITY APPROACH IN 2D 



A. Singular terms associated to plastic inclusion 

This K&M formalism can be applied to two-dimensional inclusion problems [H, [l^. Let us consider the case of a 
small inclusion of area A experiencing plastic deformation and located at the origin of the coordinate system z = 0. 
It is assumed that the stress is a constant at infinity. Outside the inclusion, the K&M potentials can be expanded as 
a Laurent scries as 



oo oo , 

^(z) = a°-z + 5: ^ , ^(z)=/3''-z + 5]% (6) 



Z" ^ Z" 



The linear terms can be easily identified as corresponding to uniform stresses while constant terms (omitted here) 
would lead to a rigid translation. It can be shown in addition that the dominant singular terms ipi/z and V'l/^ can 
be associated to the elastic stress induced by the plastic deviatoric and volumetric of an equivalent circular inclusion 
of area A. Namely, considering a circular inclusion experiencing a plastic shear strain jp and a plastic volumetric 



strain Sp we have [IE 



2i^j,Ajp 2^iASp 

n(K + 1) TT[K + 1) 



In particular, for a pure shear plastic event we obtain a quadrupolar symmetry: 

_ 2'ypfi A 



1 7rr^ 



cos(46l) (8) 



Note that we have in general to consider a complex value of 7p to include the angular dependence of the principal 
axis. In contrast, the amplitude ipi is a real number (note that the imaginary part would correspond to a point-like 
torque applied at the origin). 
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B. Generic character of the expansion 



Because of the well-known property of Eshelby circular inclusion, the above expansion limited to the ipi and 0i 
terms only is the exact (outer) solution of a uniform plastic strain distributed in the inclusion, and vanishing stress 
at infinity. However, one should note that this result is much more general. Indeed, it is seen that the physical size of 
the inclusion does not enter into the solution but through the products Ajp and A5p. Therefore, a smaller inclusion 
having a larger plastic strain may give rise to the very same field, provided the products remain constant. Therefore, 
one can consider the prolongation of the solution to a point-like inclusion (with a diverging plastic strain), as being 
equivalent to the inclusion. 

Then from the superposition property of linear elasticity, a heterogeneous distribution of plastic strain 7p(a:) in a 
compact domain, V, will give rise to such a singularity with an amplitude equal to 



and the same property would hold separately for the volumetric part. As a particular case, one finds a uniform plastic 
strain for an inclusion of arbitrary shape. 

This is the key property that allows to capture the equivalent plastic strain of an arbitrary complex configuration, 
for the above mentioned application to amorphous media. In fact, this is even the only proper way of defining the 
plastic strain for a discrete medium such as encountered in molecular dynamics simulations. The far- field behavior of 
the displacement and stress field can be accurately modelled, and without ambiguity by a continuum approach, and 
thus the above result will hold. In contrast, locally, the large scale displacement of several atoms may render difficult 
the direct computation of the equivalent plastic strain experienced in such an elementary plastic event. 

Let us however stress one difficulty: As the above argument ignores the details of the action taking place within 
the "inclusion" , plasticity has to be postulated. However, a damaged inclusion, where the clastic moduli have been 
softened by some mechanism, or even a non-linear elastic inclusion at one level of loading would behave in a similar 
way as the above plastic inclusion. Obviously, to detect the most relevant physical description, one should have 
additional informations, say about unloading. If the above amplitudes remain constant during unloading, plasticity 
would appear appropriate. If the amplitude decreases linearly with the loading, then damage is more suited. Finally, 
if the amplitudes varies reversibly with the loading, non-linear elasticity might be the best description. Thus, although 
one should be cautious in the interpretation, local damage detection from the far field may also be be tackled with 
the same tools. 



In two-dimensions, this multipole expansion formalism in the complex plane suggests to resort to contour integrals 
to extract the singularities. However the displacement field is not a holomorphic function and cannot be used directly 
for that purpose. The strategy of identification of the singularities ipn and tpn thus consists of expressing the potentials 
from the displacement field and its derivatives in order to extract the singularities via Cauchy integrals. We now simply 
express the displacement field and its derivatives: 




(9) 



C. Path independent contour integrals 
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K(p{z) - Zip'{z) - tpiz) , 



(10) 

(11) 
(12) 

(13) 



Kip'{z) - (p'{z) , 



—Zip"{z) — ij'{z) , 



This gives immediately: 
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(15) 



2/i Uj 




(16) 



Note that, except a multiplicative constant, the two last expressions are independent of materials properties. In light 
of the expansion of (/3 and ijj in Laurent scries, if an anti-clockwise contour integration is considered along a path 
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C, Cauchy residues can be formed as 



47rn(n + 1) Jc 



1 



dz 



(17) 
(18) 



These expressions can thus be obtained from the sole knowledge of the displacement field, a quantity which can be 
accessed from experiments, or from atomistic simulations. In the case of first order singularities (see Eq.([7])), the 
residue term thus only depends on the local plastic deformation (size and amplitude of deformation) and on the 
Poisson's ratio ly of the material. 

Reverting to Cartesian coordinates, where the contour is expressed as a function of the curvilinear abscissa s as 
(a;(s), y(s)), the above expression can be written 



V-i = ^j;[2[zx-y][(C/,,, 
-[i(a;2 + y2)][[/^_^^H 



~ ^V,y) ~ ^{Uy,x + Ux,y)] 



^,yy 



y,xx 



Uy.yy)]] + ^^]ds 



(19) 



IV. NUMERICAL IMPLEMENTATION 



The ultimate goal of a such a method would be to analyze numerical results obtained from molecular dynamics 
simulations of amorphous plasticity where such local structural reorganizations are expected to take place. This 
obviously raises the question of a well defined method for writing the continuous displacement field from the data of 
the discrete displacements of particles [l^] and more generally the question of the sensitivity to noise of the above 
expressions. The first point is beyond the scope of the present work and we leave it for later studies. We thus 
focus on the more restricted question of the numerical implementation and its efficiency in the case of artificially 
noise-corrupted displacement data. 

The method relies on contour integrations of derived fields of the displacement. The latter point induces a priori 
a strong sensitivity to noise. To limit such effects, first and second derivatives are extracted via an interpolation of 
the local displacement field by polynomial functions of the spatial coordinates. Moreover, the path independence of 
the contour integrals allows to perform spatial averages. We explore in the following the efficiency of this method on 
noisy data. 



A. Numerical generation of elastic fields induced by plastic inclusions 

Displacement fields are computed numerically using a finite element code, with square elements and bi-linear shape 
functions {l,x,y, xy}. Plane stress conditions of two-dimensional elasticity are used. The domain is a 150x 150 square. 
Stress free conditions are enforced all along the domain boundary. The Poisson's ratio of the material is ly = 0.20. 
Since no quantitative values of the stress are used, the value of the Young's modulus is immaterial. 

A plastic strain is implemented at the scale of one single isolated element. Within this element, the strain is the 
sum of a plastic uniform strain chosen at will, and an elastic strain. The latter is computed by solving for the two- 
dimensional elastic problem, insuring force balance and kinematic continuity at all nodes including the nodes of the 
plastic element. The chosen kinematics is too crude to solve at the scale of one single element the elastic problem 
with a good accuracy. However, remote from the inclusion, the elastic perturbation is well accounted for, and since 
all our computations are based on paths lying at a distance from the inclusion, the formalism should be applicable. A 
single element allows us to have a maximum ratio between inclusion and domain size. The price to pay for this crude 
local description is that the quantitative estimate of A'jp and A6p will differ slightly from the theoretical expectation. 
Nevertheless the path independence, (size, shape, center, ...) is expected to hold. 

We limited ourselves to such a description in order to mimic the difficulties one may face when having to deal with 
discrete element simulations. Indeed, the chosen finite element shape functions do not allow us to use this description 
strictly speaking in order to compute second order differential operators on the displacement, since the gradients of 
the latter are not continuous across element boundaries. Therefore, a regularization will be called for, as detailed 
below. 

The choice of a regular square lattice is obviously oversimplified compared with the case of the random lattices 
associated with atomistic simulations. However, as x and y directions are obviously equivalent for square elements 
and as linearity is preserved by the finite element formulation, this formulation should not introduce any breaking 
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FIG. 2: Normalized values of numerical estimates tpi and ijji obtained for three families of contours, respectively square, cross 
and rectangle shaped and of varying length L in the case of a displacement field induced by the plastic shear strain (left) or 
contraction (right) experienced by the central element of a square lattice. Theoretical expectations are ifii = i ,ipi =0 (left) 
and (fi = ,ipi — 1 (right). 



of symmetry. More specifically, the displacement field induced by a quadrupole of principal direction off axis can be 
obtained by a linear superposition of x and y components of the displacement field induced by a quadrupole aligned 
with the axis weighted by the sine and cosine of the quadrupole angle. 

Finally, the finite size of the system is also a specific difficulty encountered in practice, whereas the above argument 
uses the assumption of an infinite domain. However, such a boundary condition should not induce additional poles 
within the domain, and can be considered as a common practical difficulty encountered for all practical use of this 
tool. All these arguments are possible causes of deviation from the theoretical expectation of path independence, and 
it thus motivates a detailed study of the method stability, robustness and accuracy. 

Two test cases are studied. I: a central inclusion experiencing a shear strain 7^ = 1 (because of linearity, the actual 
amplitude is meaningless) along the x— axis; II: a central inclusion experiencing a volumetric contraction dp — —1. A 
map of the displacement fields Uy in case I is given on Fig. [T] 
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FIG. 3: Normalized values of numerical estimates ipi obtained for square contours of center (a:o,0) and of size M = 20. The 
expected behavior of 5R?/'i (unity when the inclusion lies within the contour, zero otherwise) is represented by the bold line. 
Other quantities are expected to be zero. 



B. Interpolating displacement data 

The key ingredient is to go from a continuous but non-difFerentiable displacement field obtained from the finite 
element simulation to an evaluation of the second derivative at any point in the domain. 

The results which are presented below have been obtained using the following procedure. A quadratic fit is performed 
on a square centered on one node to extract the first and second order derivatives, the obtained values are used to 
compute the integrals by quadrature. An alternative method has been tested: for an integration from (x, y) to 
(x -f 1, y), a simple fit is performed of the 12 nodes ranging from (x — 1) to (x + 2), and from (y — 1) to (y + 1), by the 
tensor product of polynomials {l,x,x'^,x^) and (l,y, y^) (12 functions). Then the integral of all required quantities 
can be computed. Estimates of derivatives using of Fourier Series with and without low pass filtering have also been 
performed. All methods give similar results provided that the area of the region used for interpolation (or filtering) 
is comparable. 



C. Path independence 

Wc first check the path independence of the contour integral in the cases I and II of isolated inclusions. For that 
purpose, we use two families of contours, square (A) and cross (B) and rectangular (C) shaped respectively as shown 
on Fig. [1] The size of these contours as well as their center can be varied. Figure [2] gives a summary of the results. 
For the three kinds of contours, we show the real and imaginary parts of the residues corresponding to Eq. (fT7|) . Note 
that the numerical results have been normalized according to the theoretical expectations ([7]) so that the expected 
numerical values are ipi = i = in case I (Fig. [2] left) and ipi = ,ipi = 1 in case II (Fig. [2] right). 

These numerical results can be considered as rather satisfactory in terms of orientation and orthogonality between 
modes ipi and ipi- the measured values of quantities whose expected value is zero remain typically below 10~^. When 
compared to their theoretical values, (^^'lear g^^^ ^contraction gxhibit relative differences of around 5-10% . Small 
fluctuations (below 5%) can be found when changing shape and size of the contours. We already commented on the 
fact that the finite element simulation are performed with a single element for the inclusion, a procedure which is 
obviously not reliable in terms of accuracy, but which allows us to have a large ratio between element and system size. 

Another test of the dependence of the numerical procedure is on the sole topology, i.e. location of the inclusion 
inside or outside the contour, we show in addition the dependence of the measured values of ipi and ipi on the location 
of the contour center. Fig. [3] shows the singularity ipi measured from the integration of displacement field II along a 
square contour centered along the x-axis. Results are normalized so that the expected value of SR'0i be unity when 
the inclusion is within the contour and zero elsewhere. The contour size is M = 20. We obtain the expected behavior: 
the values of KV'i shifts form zero to unity depending on the inclusion is within or outside the contour. Significant 
fiuctuations (10-20%) are however observed when the inclusion lies in the vicinity of the contour. 

Finally we test the dependence of the numerical method on the material properties. In the determination of cpi 
and -01 (HZl) as residues, the need to resort to a second order derivative of the displacement field is balanced by the 
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FIG. 4: Numerical estimates of \ipi\ and ipi obtained for a central inclusion experiencing a unit shear and a unit contraction 
respectively as a function of the Poisson ratio u. The numerical results obtained in plane stress conditions for a system of size 
100 X 100 with stress free boundary conditions are compared with the theoretical expectation (1 + i')/2ti. 

fact that the computation can be performed without any knowledge of the elastic properties of the material. The 
independence of the numerical procedure with the Young modulus is trivially obtained due to the linearity of the FEM 
computation. In Fig. [3]we show the dependence of the numerical results on the Poisson ratio. FEM computations have 
been performed on systems of size 100 x 100 with stress free boundary conditions and a central inclusion experiencing 
a unit shear and a unit contraction respectively. Poisson ratios have been varied from 0.05 to 0.45 by step of 0.05. 
The results shown in the figure compare the numerical estimates obtained for a square contour of size 40 centered on 
the inclusion with the theoretical expectation ipi ^ Lpi = {\ + j/)/27r. The numerical results show that the volumetric 
strain is weakly dependent on Poisson's ratio, but the elementary shear is more poorly estimated for high Poisson's 
ratio. 



V. DISCUSSION/CONCLUSION 



The proposed approach is based on a an exact result and hence theoretically establish a parallel with other types 
of elastic singularities (in particular for stress intensity factors which characterize crack loadings) where similar path 
integral are well known. When tested on direct numerical simulations, we could recover the main topological properties 
expected in this context: path independence and detection of the absence/presence of a singularity within the contour. 
However, the quantitative results proved more disappointing: the method is rather unprecise on the determination 
of the prefactor of the singularity and is more generally rather sensitive to noise. The main cause is presumably 
due to the inconsistent regularity of the displacement field solution (simple continuity) with the need to resort to 
estimates of first and second order differentials. A piecewise high order polynomial interpolation is operational for 
integrals over finite segments, however, from one segment to the next, first and second order differentials will display 
a discontinuous character, which obviously affect the method and result. Moreover, being a path integral, the method 
does not take advantage of the knowledge of the displacement field at all points of a domain. To make the method 
more robust with respect to noise, different approaches can be followed. One natural way is to average the result over 
different contours, thus transforming the contour integral into a domain integral. An arbitrary weight average can 
also be considered, and hence, one could optimize the weight in order to achieve the least noise sensitivity. Such a 
method was explored successfully for cracks in ref. [l3|. Note finally that extensions to three dimensions obviously 
require a different technique than Kolossov and Muskhelishvili potentials, and contour integrals, however still a linear 
extraction operator acting on the displacement field can be computed to provide similarly the equivalent plastic strain. 
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